#include <iostream>
#include <opencv2/opencv.hpp>
#include <Eigen/Dense>
#include <opencv2/core/eigen.hpp>

// 创建高斯核
Eigen::MatrixXf createGaussianKernel(int kernelSize, float sigma) {
    Eigen::MatrixXf kernel(kernelSize, kernelSize);
    float sum = 0.0;
    int halfSize = kernelSize / 2;

    for (int i = -halfSize; i <= halfSize; ++i) {
        for (int j = -halfSize; j <= halfSize; ++j) {
            float value = std::exp(-(i * i + j * j) / (2 * sigma * sigma));
            kernel(i + halfSize, j + halfSize) = value;
            sum += value;
        }
    }
    kernel /= sum;
    return kernel;
}

// 计算频域特征
void calculateFrequencyFeatures(const cv::Mat& magnitudeImg, double& lowFreqEnergy, double& highFreqEnergy) {
    int cx = magnitudeImg.cols / 2;
    int cy = magnitudeImg.rows / 2;
    int radius = std::min(cx, cy) / 2;

    lowFreqEnergy = 0.0;
    highFreqEnergy = 0.0;

    // 确保 magnitudeImg 的数据类型是 CV_32F
    CV_Assert(magnitudeImg.type() == CV_32F);

    for (int i = 0; i < magnitudeImg.rows; ++i) {
        for (int j = 0; j < magnitudeImg.cols; ++j) {
            float value = magnitudeImg.at<float>(i, j); // 使用 float 类型访问
            int dist = std::sqrt((i - cy) * (i - cy) + (j - cx) * (j - cx));
            if (dist < radius)
                lowFreqEnergy += value;
            else
                highFreqEnergy += value;
        }
    }
}

// 可视化 DFT
cv::Mat visualizeDFT(const cv::Mat& img, cv::Mat& magnitudeFloat) {
    // 将输入图像转换为复数形式
    cv::Mat planes[] = {cv::Mat_<float>(img), cv::Mat::zeros(img.size(), CV_32F)};
    cv::Mat complexImg;
    cv::merge(planes, 2, complexImg);

    // 计算 DFT
    cv::dft(complexImg, complexImg);

    // 分离实部和虚部
    cv::split(complexImg, planes);
    cv::magnitude(planes[0], planes[1], magnitudeFloat);

    // 对幅值取对数并归一化
    cv::log(magnitudeFloat + 1, magnitudeFloat);
    cv::normalize(magnitudeFloat, magnitudeFloat, 0, 1, cv::NORM_MINMAX);

    // 创建一个副本用于伪彩色处理
    cv::Mat magnitude8U;
    magnitudeFloat.convertTo(magnitude8U, CV_8U, 255);

    // 应用伪彩色映射
    cv::Mat magnitudeColorMap;
    cv::applyColorMap(magnitude8U, magnitudeColorMap, cv::COLORMAP_JET);

    return magnitudeColorMap;  // 返回伪彩色图像
}


int main() {
    // 加载图片
    cv::Mat img = cv::imread("../pic.png", cv::IMREAD_GRAYSCALE);
    if (img.empty()) {
        std::cerr << "无法加载图片！" << std::endl;
        return -1;
    }
    cv::resize(img, img, cv::Size(512, 512));

    // 创建高斯核
    Eigen::MatrixXf gaussianKernel = createGaussianKernel(5, 10.0);
    cv::Mat gaussianKernelCV;
    cv::eigen2cv(gaussianKernel, gaussianKernelCV);

    // 应用高斯模糊
    cv::Mat blurredImg;
    cv::filter2D(img, blurredImg, -1, gaussianKernelCV);

    // 显示原始图像与模糊图像
    cv::imshow("Original Image", img);
    cv::imshow("Blurred Image", blurredImg);

    // 计算频域特征
    double originalLowFreq, originalHighFreq;
    double blurredLowFreq, blurredHighFreq;

    cv::Mat originalSpectrumFloat, blurredSpectrumFloat;
    cv::Mat originalSpectrum = visualizeDFT(img, originalSpectrumFloat);
    cv::Mat blurredSpectrum = visualizeDFT(blurredImg, blurredSpectrumFloat);

    calculateFrequencyFeatures(originalSpectrumFloat, originalLowFreq, originalHighFreq);
    calculateFrequencyFeatures(blurredSpectrumFloat, blurredLowFreq, blurredHighFreq);
    
    // 显示频域图像
    cv::imshow("Original Spectrum", originalSpectrum);
    cv::imshow("Blurred Spectrum", blurredSpectrum);

    // 输出频域特征
    std::cout << "Original Image Frequency Characteristics:\n";
    std::cout << "  Low Frequency Energy: " << originalLowFreq << std::endl;
    std::cout << "  High Frequency Energy: " << originalHighFreq << std::endl;
    std::cout << "  Low-to-High Frequency Ratio: " << originalLowFreq / originalHighFreq << std::endl;

    std::cout << "\nBlurred Image Frequency Characteristics:\n";
    std::cout << "  Low Frequency Energy: " << blurredLowFreq << std::endl;
    std::cout << "  High Frequency Energy: " << blurredHighFreq << std::endl;
    std::cout << "  Low-to-High Frequency Ratio: " << blurredLowFreq / blurredHighFreq << std::endl;

    cv::waitKey(0);
    return 0;
}